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Wang and Landau proposed recently, a simple and flexible non-Boltzmann Monte 
Carlo method for estimating the density of states, from which the macroscopic prop- 
erties of a closed system can be calculated. They demonstrated their algorithm 
by considering systems with discrete energy spectrum. We find that the Wang- 
Landau algorithm does not perform well when the system has continuous energy 
spectrum. We propose in this paper modifications to the algorithm and demonstrate 
their performance on a lattice model of liquid crystalline system (with Lebwohl- 
Lasher interaction having continuously varying energy), exhibiting transition from 
high temperature isotropic to low temperature nematic phase. 

PACS numbers: 05.10.Ln; 61.30.pq; 64.70.Md 

I. INTRODUCTION 

Monte Carlo methods have emerged as a powerful and reliable tool for simulating several 
complex phenomena in statistical physics, see e.g. Q, Q|- The Metropolis algorithm [3] 
discovered in the middle of the last century can be considered as the starting point. This 
algorithm generates a Markov chain, the asymptotic part of which contains microstates 
belonging to a canonical ensemble at a temperature chosen for the simulation. Expectation 
value of a macroscopic property can be estimated by taking a simple arithmetic average 
over a Monte Carlo sample. The associated statistical error is inversely proportional to the 
square root of the sample size. Thus, in principle, we can estimate a physical property to the 
desired accuracy by simply increasing the sample size. However, if successive microstates in 
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the sampled Markov chain are correlated, the statistical error increases by a factor a/1 + 2r*, 
see e.g. 0], where r* is the integrated correlation time. Such a situation obtains when we 
simulate a system close to criticality. 

The Metropolis algorithm and its several variants like Glauber heat-bath [6] and 
Kawasaki exchange pj algorithms come under the class of Boltzmann sampling techniques. 
The limitations of Boltzmann sampling have long since been recognized. For example it 
can not address satisfactorily problems of critical slowing down, i.e. divergence of r* with 
increase of system size, near continuous phase transition. Cluster algorithms ||| overcome 
this problem. Boltzmann sampling is also not suitable for problems of super critical slowing 
down near first order phase transitions. The microstates representing the interface between 
ordered and disordered phases have intrinsically low probability of occurrence in a closed 
system and hence are scarcely sampled; switching from one phase to the other takes a very 
long time due to the presence of high energy barriers when the system size is large; as a result 
the relative free energies of ordered and disordered phases can not be easily and accurately 
determined. Finally it is quite difficult to estimate the absolute values of entropy or free 
energies in Boltzmann sampling techniques. 



A. Non-Boltzmann sampling 



That non-Boltzmann sampling can provide a legitimate and often superior alternative to 
Boltzmann sampling was recognized even during the early days of Monte Carlo practice, see 
e.g. [9]. However, the practical convenience and significance of non-Boltzmann sampling 
was appreciated only in the middle of seventies when Torrie and Valleau proposed the so 
called umbrella sampling; this is a fore-runner to all the subsequent non-Boltzmann sampling 
techniques includin g th e multicanonical Monte Carlo and its several and recent variants. 
Entropic sampling , equivalent to multicanonical sampling provides a transparent 
and intuitively appealing insight into non-Boltzmann Monte Carlo techniques. It is based 
on the following premise. 

The probability that a closed system can be found in a microstate C is given by 



P(C) = Z(f3) 



-i 



exp 



PE(C) 



(1) 



In the above E(C) is the energy of the microstate C and (3 = (ksT) 1 , where ks is the 
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Boltzmann constant and T is the temperature. Z(f3) is the canonical partition function 
given by, 



Z{(3) = £exp[-/?£(C) 



D(E)exp{-/3E)dE, 

where D(E) is the density of states. The probability density for a closed system to have an 
energy E is given by, 



P B {E) oc D(E)exp(-pE) 



(2) 



where the suffix B indicates that it is the Boltzmann-Gibbs distribution (appropriate for 
modeling a closed system). Let us suppose we want to sample microstates in such a way 
that the resultant probability density of energy is given by, 

-l 



P g (E) ex D(E) g(E) 



(3) 



where g(E) is chosen as per the desired non-Boltzmann distribution. Once P g {E) is known, 
an ensemble consistent with Eq. (JHJ) can be constructed as follows. 

Let Ci be the current microstate and Ct the trial microstate obtained from Ci by making 
a local change, e.g. flip a randomly chosen spin in Ising model simulation. Let E^ = E(Ci) 
and E t = E(C t ) denote the energy of the current and of the trial microstate, respectively. 
The next entry Q+i in the Markov chain is taken as, 



C 



C t with probability p, 



i+l 



(4) 



Ci with probability (1 — p), 
where the acceptance probability p is given by, 



p = mm 





= min 


L g(Ei)] 


[ ' P 9 {Ei). 




[ ' g{Et)\ 



(5) 



We call this non-Boltzmann sampling with respect to a given g{E). It is easily verified 
that the above acceptance rule obeys detailed balance. Hence the Markov chain constructed 
would converge asymptotically to the desired ^-ensemble. 



4 



When [^(i?)] -1 = exp(— (3E) we recover conventional Boltzmann sampling, implemented 
in the Metropolis algorithm. For any other choice of g(E) we get the corresponding non- 
Boltzmann sampling. Now, canonical ensemble average of a macroscopic property 0(C) 
can be obtained by un- weighting and re- weighting of 0(C) for each C sampled from the 
^-ensemble. For un- weighting we divide by [^(^(C))] -1 and for re- weighting we multiply by 
exp[—(3E(C)]. The weight factor associated with a microstate C belonging to the ^-ensemble 
is thus, 



W(C, 0) = g (JS?(C)) exp [ - (3E(C)] . (6) 



We then have, 



E c 0(C)^(C,/3) 

The left hand side of the above is the equilibrium value of O in a closed system at /3, while on 
the right side the summation in the numerator and in the denominator runs over microstates 
belonging to the non-Boltzmann g-ensemble. It is also clear that from a single simulation 
of a g-ensemble, we can calculate the canonical average of O at various temperatures. 



B. Entropic sampling 

Entropic sampling obtains when g(E) = D(E). This choice of g(E) renders P g {E) the 
same for all E, see Eq. The system does a simple random walk on a one dimensional 
energy space. Hence all energy regions are visited with equal probability. As a result, 
in the case of first order phase transition for example, the microstates on the paths (in 
the configurational space) that connect ordered and disordered phases would get equally 
sampled. A crucial issue that remains to be clarified pertains to the observation that we do 
not know D(E) a priori. 

In entropic sampling we employ a strategy to push g(E) closer and closer to D(E), 
iteratively. We divide the range of energy into a large number of bins of equal widths. We 
denote the discrete-energy version of g(E) by the symbol {#j : % = 1, 2, • • • }. We start with 
{gf 1 ' = 1 V i}; the superscript is iteration run index and the subscript is energy bin index. 
The aim is to update {^} from one iteration to the next: {gf*} — > {g^} —>•••• {g[ } 
so that asymptotically we get {gi\ as close to {A} as desired, where {A} is the discrete 
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energy representation of D(E). The iteration is carried out as follows. In the k — th 
iteration, for example, we generate a large number of microstates employing acceptance 
probability based on {<?| } and accumulate an histogram {hi : i — 1, 9L 3 } of energy of 
visited microstates. We update {g\ : i — 1, 9L 3 } to {gf^ : i — 1, 9L 3 }, as given below, 

' gf ] if^ = o, 

gt 1] = (8) 
x hi if hi ^ 0, 

for all i = 1, 2, ■ • ■ , 9L 3 . The updated {gf - ^} is employed in the next i.e. (k + 1) — i/i 
run, during which a fresh histogram of energy is generated. After each run, the histogram 
is examined for its uniformity. Flatter the histogram, closer is {g{\ to {-D«}. Thus, the 
calculated histogram serves two purposes in entropic sampling, one for updating and 
the other for monitoring the convergence of {g{\ to {D{\. However, it is often neither 
practical nor necessary to get a strictly flat histogram; an approximately flat histogram 
would be adequate, thanks to the un-weighting followed by reweighting with the Boltzmann 
rule while calculating the averages, see Eqs. (0 0). Hence the calculated macroscopic 
properties would come out right, even if {g{\ does not converge strictly to {-Dj}. 



C. Wang-Landau algorithm 

A simple and flexible variant to entropic sampling was proposed recently by Wang and 
Landau [13]. The distinguishing feature of this algorithm is the dynamic evolution of the 
acceptance probability, p; we update {^} after every Monte Carlo step. Let us say the 
system visits a microstate in a Monte Carlo step and let the energy of the visited microstate 
fall in the m-th energy bin; then g m is updated to / x g m , where / is the Wang-Landau 
factor, see below. The updated {gi} becomes operative immediately for determining the 
acceptance/rejection criteria of the very next trial microstate. We set / = / for the zeroth 
run; fo can be any number greater than unity; the choice of fo = e has been originally 
recommended by Wang and Landau. We generate a large number of microstates employing 
the dynamically evolving p. At the end of a run we calculate the histogram of energy of 
microstates visited by the system during the run. Because of the continuous updating of 
p, the energy span of the density of states increases significantly and the energy histogram 
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serves to monitor the convergence of {<&} to {-Dj}. A run should be long enough to facilitate 
the system to span the energy over the desired range and to render the histogram of energy 
approximately flat. At the end of, say, the u-th run, the Wang-Landau factor for the next 
run is set as / = f v+ i = \pf v - After several runs, this factor would be very close to unity; this 
implies that there would occur no significant change in {g{\ during later runs. For example 
with the square-root rule and f = e, we have / 25 = exp(2~ 25 ) < 1 + 10~ 7 . It is clear that 
/ decreases monotonically with increase of the run index and reaches unity asymptotically. 
Wang and Landau have recommended the square- root rule 13|; any other rule consistent 
with the above properties of monotonicity and asymptotic convergence to unity should do 
equally well. 

From the converged g the desired macroscopic properties of the system can be calculated; 
to this end we invoke the the connection between the density of states and microcanonical 
entropy, a(E) = ks\ogD(E). Thus the Monte Carlo estimate of microcanonical entropy is 
&b logg(E'). For implementing such a scheme we need to normalize g{E). The normalization 
constant should be obtained from known properties of the system. For example in Ising 
model, the ground state is doubly degenerate: D(E m[n ) = 2. The total number of microstates 
equals 2 V where V is the number of Ising spins in the Monte Carlo model: fjf max D(E)dE = 

^ -^min 

2 V . Either of these known information can be employed for normalizing g. The normalized 
g(E) provides a good approximation to D(E). 

Alternately, we can take the output {g{\ from the above and carry out a single long non- 
Boltzmann sampling run which generates microstates belonging to the g-ensemble. (Note 
that during the production run we do not update g(E). By un-weighting and re-weighting 
of the microstates generated in the production run, we calculate the desired properties of 
the system as a function of (3. This is the strategy we shall follow for the simulation of liquid 
crystal system, described in the rest of the paper. In this strategy, we can employ arbitrary 
normalization of g; more importantly, it is adequate if 

(a) the system visits the energy region of interest and not necessarily the entire range and 

(b) the histogram of energy in the region of interest is approximately flat. 

The usefulness of the Wang-Landau algorithm has been unambiguously demonstrated for 
systems with discrete energy spectrum. However, when we try to apply this technique to 
systems with continuous energy, there are serious difficulties. Liquid crystalline materials 
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with continuous energy spectrum provide such an example. In this paper we report simu- 
lation of a liquid crystalline system focusing attention on nematic-isotropic transition. The 
paper is organized as follows. In section II, we describe a lattice model and Hamiltonian 
of a liquid crystal system. We simulated the system with the conventional Wang-Landau 
algorithm. We find that the dynamics becomes extremely slow even for moderately large 
systems. The system gets stuck in certain regions of the configurational space. This problem 
appears to be generic to the algorithm when applied to continuous energy systems. Hence 
we modify the Wang-Landau algorithm and the details of the simulation are given in section 
III. The results on temperature variation of various macroscopic properties of the system 
are discussed in section IV. In section V we briefly summarize the work and highlight its 
salient features. 



II. LATTICE MODEL OF BULK LIQUID CRYSTALS 

We consider an L x L x L cubic lattice with each lattice site holding a three dimensional 
unit vector \u), called a spin. The elements of the vector are the direction cosines of a 'spin' 
in a laboratory frame of reference. A 'spin' represents notionally, a single uniaxial liquid 
crystal molecule or more realistically a cluster containing typically a hundred of them. The 
spins are actually 'headless' in the sense the system has head-tail flip symmetry. Two nearest 
neighbour spins interact with each other as per a potential proposed by Lebwohl and Lasher 
(LL) 14] which has such a head-tail flip symmetry. The interaction energy is given by, 

e^ = -^[3cos 2 (^)-l], (9) 

where i and j are nearest neighbour lattice sites. 9^ is the angle between the two spins: 
cos(8ij) = (ui\uj). The interaction energy of a single nearest neighbour pair of spins ranges 
from a minimum of —1, when 9{j = 0, or equivalently (v,i\uj) = 1, to a maximum of +1/2, 
when 9ij = tt/2 or equivalently (ui\uj) = 0. Total energy of the system in microstate C is 
given by, 

E(C) = J2%v (10) 

M 

where the sum runs over all distinct nearest neighbour pairs of lattice sites in the system 
taking into account the periodic boundary conditions in all the three directions. The total 
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energy of the system thus varies continuously from a minimum of — 3L 3 to a maximum of 
+3L 3 /2. When the system is completely ordered with all the spins aligned, the energy is 
minimum and equals — 3L 3 ; the energy is zero for an isotropic (completely disordered) phase. 
We calculate several macroscopic properties of the liquid crystalline system and report here 
results which include orientational order parameter (S), average energy (E), specific heat 
(at constant volume) Cy, and the Binder's reduced fourth cumulant of energy V4. 

First we employed conventional Wang-Landau algorithm and carried out Monte Carlo 
simulation of the lattice model of the Liquid crystalline system. We found that the dynamics 
was extremely slow when the system size L is 6 and above. The calculated density of states 
g(E) becomes steeper with increase of Monte Carlo iterations. As a result the system gets 
stuck in certain narrow regions. There is practically no evolution of the calculated density of 
states g(E). Increasing the number of Monte Carlo steps in a Wang-Landau iteration does 
not seem to remedy the situation. Instead, sharp peaks emerge and grow at either ends of 
g(E). These problems appear to be generic to the Wang-Landau algorithm when applied to 
continuous energy systems. To overcome them, we experimented with several modifications 
[lfij ] of the algorithm and finally arrived at a strategy described in the next section. 



III. MODIFIED WANG-LANDAU MONTE CARLO SIMULATION OF BULK 

LIQUID CRYSTAL SYSTEM 

Free wheeling spins are placed on the vertices of a three dimensional cubic lattice with 
their orientations sampled randomly and independently. Orientation of a spin is specified 
by the polar angle 6 and an azimuthal angle (f) with respect to a laboratory fixed three 
dimensional co-ordinate system. We sample fi = cos(#) from a uniform distribution between 
and 1, and the azimuthal angle from a uniform distribution between and 2n. We divide 
the energy range (— 3L 3 , +1.5L 3 ) into 9L 3 number bins of equal widths AE = 0.5. We start 
with an array {gi = e 2 V % = 1, 9L 3 }. Let Co denote the initial microstate constructed by 
placing spins with random and independent orientations at the lattice sites. Let E(C ) fall 
in the [i — th energy bin. We select randomly a lattice site and make a random change in the 
orientation of the spin residing at that site. We employ Barker's method ^3] to generate 
a new trial orientation from the current microstate. Let Ct denote the trial microstate and 
let its energy belong to the ^-th bin. If g v < g^, we accept the trial state and set C^ = C t ; 
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if g v > g^, we calculate the ratio p = g^,/g u - We select a random number r (uniformly 
distributed between to 1); if r < p we accept the trial microstate: C^ = C t . Otherwise we 
reject Cf and set C[ = Co- This constitutes a single move. We continue in the same fashion 
and get, Co — > — > C 2 • • • C' L3l — > C±. A set of L 3 moves constitutes a Monte Carlo Sweep 
(MCS). In the first MCS, since gi is the same for all i, every move gets accepted. 

At the end of the MCS, let us say the system is in microstate Ci. Let E{C\) belong to 
the k-th energy bin. We update <j% to / x g k , where / = f . The updated {gi} becomes 
operative for deciding acceptance/rejection in the next L 3 moves that constitute the next 
MCS leading to C2. Thus we get a chain of microstates Co — > Ci — > •••Cat. We take 
N = 10, 000. Generating a Markov chain of length iV constitutes one iteration. For the next 
iteration we change / to / 9 . The microstate generated at the end of an iteration forms the 
initial microstate for the next. Also the updating of the density of states is continued from 
one iteration to the next. We carry a total of M iterations and this constitutes a Wang- 
landau (W-L) run, see below. In the last iteration of a W-L run, we have / = Jm = fo , 
where /i(M) = (0.9) M . We have chosen M = 160 so that f M -lm 10~ 7 for /„ = 10. 

We start a W-L run with / reset to fo. The microstate generated at the end of a W-L 
run is taken as the initial microstate for the next. Similarly the updated density of states 
{<7j} at the end of a W-L run provides the initial density of states for the next. We carry 
out a total of 50 W-L runs. The value of fo is 100 for the first forty, 10 for the next 9 and 
fo = e for the last W-L run. 

The density of states at the end of 50 W-L runs is taken as an input for a long non- 
Boltzmann sampling run of 2.5 million Monte Carlo sweeps, called the production run. 
Thus we get a g— ensemble of microstates from which the desired macroscopic properties 
can be calculated by un-weighting and re-weighting. 

We also found that it is imperative to employ numerical techniques that avoid overflow 
problems and the attendant loss of precision due to truncation. To this end, we adapted the 
techniques suggested by Berg [l^. These involve principally the following. Let aij = log^j 
denote the microcanonical entropy. We define & = loga«. We carry out all the calculations 
in terms of {£j : i — 1, 2, • • ■ , 9L 3 }. We derive expressions for acceptance probability p 
in terms of and employ them in the simulation. Similarly we derive expressions for 
the updating of and for un-weighting and re-weighting, in terms of These are 

briefly described in the appendix. Employing this modified Wang-Landau algorithm we 
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simulated a lattice model of liquid crystalline system with L = 4, 6, 8, 10 and 12 focusing on 
nematic-isotropic transition. We present the results in the next section. 

IV. CHARACTERISTICS OF BULK LIQUID CRYSTALS 

The orientational order parameter S of a microstate is defined as follows. Let \ui) denote 
the spin at the 2-th lattice site. We first construct an average projection operator for a 
microstate C given by, 

L 3 

^(C) = ^f>*>(^l ■ ( H ) 

i=l 

From A we construct a traceless symmetric tensor, 

Q(C) — A — itrace(A) x I , (12) 

where / denotes a unit matrix. Let X max (C) denote the largest eigenvalue of Q. Then 
S(C) = 3\ max (C)/2. The corresponding eigenvector |A max ) defines the director for the 
microstate C. 

Figure (JTJ) depicts (S) as a function of temperature for system sizes L = 4, 6, 8, 10 and 
12. We observe that the modified Wang-Landau Monte Carlo simulation predicts correctly 
the transition from a high temperature disordered (isotropic) phase to a low temperature 
nematic phase. For L = 4 the transition is not sharp, due to finite size effects. However 
when we increase the system size, the transition becomes sharper. 

Next we investigate the behaviour of specific heat at constant volume Cy as a function 
of temperature. Cy is calculated from energy fluctuations and is given by, 

C - {E2) ~ {E)2 (13) 
° V ~ k B T* ■ (13) 

The results are depicted in Fig. (J2J). As L increases the Cy profile becomes sharper. Also the 
temperature T^i(L) at which the specific heat is maximum, shifts slightly to lower values, 
as expected. The transition temperature for the LL model has been obtained earlier [19j | 
and is given by T^i = 1.1237 ± 0.0006. We find that Tni(L = 12) calculated from the C r V 
is 1.126, in good agreement with the earlier estimate. We can estimate the L — > oo limit of 
the transition temperature by finite size scaling discussed later. 
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The variation of average energy (E) with temperature is depicted in Fig. (jSJ). This 
quantity decreases with decrease of temperature. At transition the fall is sharp for large L. 

We have calculated Binder's reduced fourth order cumulant of energy denoted by the 
symbol V4, see It is given by, 

n = i- JFl. (14) 

3(E 2 ) 2 

The variation of V4 with T is depicted in Fig. (0J) for L = 4, 6, 8, 10 and 12. Each 
curve shows a minimum at an effective transition temperature. As the system size increases 
the effective transition temperature shifts to lower values as expected. Also V4 — > 2/3 for 
T << Tni for all L considered and for T » Tni for large L. This is a a clear signature of 
a first order transition. 

Figure (jSJ) depicts microcanonical entropy a(= logg) versus E for system size L = 12 on 
a log-linear graph. The curves depict the shape of £(loga) versus E. The results on entropy 
after successive W-L runs are shown starting from the inner most and ending in the outer 
most. We see clearly that the range of energy spanned increases with increase of W-L runs. 
The outer most curve is the output of the last W-L run.. The data on : i = 1,9L 3 } 
obtained at the end of the last W-L run is employed in the long non-Boltzmann sampling run 
(production run) and a g-ensemble of microstates is generated. All the quantities referred 
to above were calculated by un-weighting and re-weighting at temperatures spaced out with 
a fine resolution of 0.001. 

Finally we have presented in Fig. © the finite size scaling of the transition temperature 
obtained from specific heat, orientational susceptibility and the fourth order cumulant of 
Binder. The orientational susceptibility \ is calculated from the fluctuations of 5* and is 
given by, 

X ~ k B T ■ (15) 

The transition temperature is plotted against inverse of the volume of the system. The 
three curves scale linearly with 1/L 3 and extrapolation [L — > 00) gives an estimate of the 
nematic-isotropic transition temperature. T^i(L = 00) estimated from specific heat data is 
1.1284, from susceptibility data is 1.1299 and from the data on Binder's cumulant is 1.1211. 
These results are in good agreement with Tni = 1.126(5) - an earlier estimate, see |3] upto 
second decimal. 
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V. CONCLUSIONS 



We have demonstrated for the first time the applicability of the recently proposed Wang- 
Landau Monte Carlo algorithm to the study of liquid crystalline systems with continuous 
energy spectrum. We have made use of the flexibility of the algorithm and studied nematic- 
isotropic transition in a three dimensional lattice model of bulk liquid crystals. We have 
employed the Lebwohl-Lasher potential that has the head-tail flip symmetry of the nematic 
director and in which the energy varies continuously. For even moderately large system 
Wang-Landau dynamics becomes unacceptably slow. The density of states function g(E) 
gets confined to a narrow energy range and becomes steep. As a result, the system spans 
only a restricted range of energy. Increasing the number of sweeps in an iteration does not 
seem to help. This slowing down of dynamics seems to be an inherent problem of this algo- 
rithm for such systems. Interestingly such problems do not arise for simulating systems with 
discrete energy spectrum e.g. Ising and Potts spin models. To overcome these problems, 
we have proposed a few modifications to the Wang-Landau algorithm and demonstrated 
that these modifications help the basic algorithm to span larger regions of the energy space 
systematically. We show that the macroscopic properties bulk liquid crystalline system can 
be calculated with a good degree of accuracy and with vastly improved temperature reso- 
lution. This opens up the possibility of exploiting the full power of the (non-Boltzmann) 
Wang-Landau Monte Carlo techniques to simulate several complex phenomena in liquid 
crystalline systems. Examples of such problems include phase transition in thin films de- 
posited on substrates with complex geometry, study of polymer dispersed liquid crystals 
(PDLC) and effect of disorder and confinement on the nematic-isotropic phase transition. 
Work on these and related problems are in progress and will be reported soon. 
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Appendix 



Let g>j denote the number of microstates in the i-th energy bin and oti = log(<7j) the 
corresponding microcanonical entropy. We define & = log{oti). In the program only the 
array {£, : i = 1,9L 3 } is stored, updated and eventually employed in reweighting. All the 
required parameters like the acceptance probability p, and un-weighting cum re-weighting 
factor W are calculated in terms of : i — 1, 9L 3 }. 

First we initialize {£ = log(2) V % — 1,9L 3 }. Let the energy of the current microstate 
belong to an energy bin c and the trial microstate, to an energy bin t. The acceptance 
probability of the trial state in the Wang-Landau algorithm is given by by, 



V 



1 if it < £ c 



cxp 



- exp | & + log (l - exp (- (& - £.)) 



if ic < & 



If the visited microstate has an energy falling in the say i-th bin, then £j is updated to 
ii + log(log(/)) where / is the Wang-Landau factor for that run. 

The un-weighting and re- weighting of microstates belonging to the g— ensemble is carried 
out as follows. Let C be the microstate under consideration. Let the energy E = E(C) of 
the microstate fall in the bin c. Let £ c be the value of £ in that bin. The weight factor W(C) 
attached to C G g — ensemble is given by 



W{C) 



exp 



+ exp 



cxp 



— exp 



£ c + log{l-exp(-A 1 )} 



log(/?£0 + log{l-exp(-A 2 )} 



where A x = £ c - log(/3£) > 



where A 2 = log(fiE) - £ c > 



It is easily verified that the above weight factor is identical to the one given by Eq. (JHJ) except 
that we have expressed it in terms of instead of {gi}- The average of a macroscopic 
property O is calculated employing Eq. (J2J). 
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FIG. 1: The orientational order parameter S versus temperature for L = 4, 6, 8, 10 and 12. The 
transition becomes sharper with increase of system size 




FIG. 2: Specific Heat Cy as a function of temperature for L = 4, 6, 8, 10, andl2; the transition 
becomes sharper with increase of system size; the transition temperature (the value of T at which 
the curve peaks) shifts to lower values with increase of system size. 




FIG. 3: Average energy as a function of temperature for L = 4, 6, 8, 10, andl2, from bottom to 
top respectively. The transition becomes sharper with increase of system size 




FIG. 4: Binder's fourth order cumulant of energy for L = 4, 6, 8, 10, and 12. The behaviour is 
indicative of first order transition 
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FIG. 5: Microcanonical entropy ccj = logt^ as a function of energy at the end of successive outer 
iterations starting from the inner most curve to the outer most. The data correspond to L = 12 
and are plotted on a log-linear curve. The shape of the curve will correspond to that of x employed 
in the simulation program. The logarithm of the outermost curve is taken as the input for a 
long production run which generates a g— ensemble. Note that almost the entire energy range is 
spanned. 
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FIG. 6: The transition temperature versus 1/L 3 . The top line with symbol A corresponds to 
Tni obtained from fourth cumulant of Binder; the middle line with symbol □ corresponds to Tni 
obtained from the orientational susceptibility and the bottom line corresponds to the Tnj obtained 
from the specific heat. The value of Tni in the limit L —> oo is given by 1.1284 from finite size 
scaling of specific heat data, 1.1299 from data on susceptibility and 1.1211 from data on Binder's 
fourth cumulant of energy. 



